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We derive expressions for the lowest nonlinear elastic constants of amorphous solids in ather- 
mal conditions (up to third order), in terms of the interaction potential between the constituent 
particles. The effect of these constants cannot be disregarded when amorphous solids undergo in- 
stabilities like plastic flow or fracture in the athermal limit; in such situations the elastic response 
increases enormously, bringing the system much beyond the linear regime. We demonstrate that the 
existing theory of thermal nonlinear elastic constants converges to our expressions in the limit of 
zero temperature. We motivate the calculation by discussing two examples in which these nonlinear 
elastic constants play a crucial role in the context of elasto- plasticity of amorphous solids. The first 
example is the plasticity-induced memory that is typical to amorphous solids (giving rise to the 
Bauschinger effect). The second example is how to predict the next plastic event from knowledge 
of the nonlinear elastic constants. Using the results of this paper we derive a simple differential 
equation for the lowest eigenvalue of the Hessian matrix in the external strain near mechanical in- 
stabilities; this equation predicts how the eigenvalue vanishes at the mechanical instability and the 
value of the strain where the mechanical instability takes place. 



I. INTRODUCTION 

Many phenomena occurring in elastic materials, 
like plasticity, fracture and shear banding were tradi- 
tionally studied assuming that the framevifork of linear 
elasticity can be employed to describe the dynamics of 
the evolving systems. Recently it became clearer that 
this is not a good idea; close to the fracture tip, where 
the stress field tends to diverge, or near a plastic in- 
stability, where the shear modulus diverges towards 
—00, nonlinear effects become crucial if not dominant 
Moreover, recent studies indicate that non- 
linear elastic moduli play an important role in the 
memory that amorphous solids exhibit of their load- 
ing trajectory Q . The well known Bauschinger effect 
can be explained as a result of the growth of the sec- 
ond order elastic modulus which is identically zero 
in an isotropic amorphous solid. It becomes there- 
fore necessary to present a microscopic theory of the 
nonlinear elastic constants to the aim of computing 
them in numerical simulations. The present paper 
has in mind athermal quasistatic simulations, a sub- 
ject of high theoretical interest for which there had 
been great recent progress in understanding the range 
of phenomena observed, including a variety of scaling 
laws and the emergence of subsequent scaling theories 
[^-Q. Interestingly, the finite-temperature counter- 
part of the theory presented in this paper is available 
in the literature, but it is not quite obvious how to 
extract from it the athermal limit. We will show be- 
low that indeed the thermal theory converges to our 
theory when T — )■ 0. 

The structure of this paper is as follows. In Sec. |ll] 
we present the derivation of the athermal elastic con- 
stants. In Sec, mil we present a derivation of the ther- 
mal elastic constants; the results of this derivation 
are scattered in the literature, but it is worthwhile to 
present them here in a compact and consistent nota- 



tion. In Sec. II VI we demonstrate that the limit T — >■ 
of the thermal elastic constants reduces properly to 
our results in Sec. [TTl Finally, in Sec. |V] we demon- 
strate the usefulness of our results in the context of 
the plasticity induced memory of amorphous solids 
and in predicting mechanical instabilities. 



II. DERIVATION OF THE ATHERMAL 
NONLINEAR ELASTIC CONSTANTS 

A. definitions 

We denote the v component of the position of the 
i'th particle as x]^\ in the following Roman characters 
denote particle indices, and Greek characters denote 
cartesian components. Given a linear transformation 
of coordinates ffi/^xj,, the resulting displace- 

ment field is u\, — H^f^xl. — x]^. The strain tensor is 
defined up to second order in the derivatives of the 
displacement field as 
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where here and below repeated indices are summed 
over, unless explicitly indicated otherwise. In terms 
of the transformation H, the strain tensor is 



Eaf! = \ {Hg^Hi^p - Safs) . (2) 

The displacement field in the definitions ([T]) and ([2]) 
should be understood as taken with respect to the ac- 
tual configuration, which may be arbitrarily deformed; 
the identity of the un-deformed isotropic reference 
state is of no interest in the following derivation. 

The free energy T is defined in terms of the parti- 
tion function Z as J- = — TlogZ, where 
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and E = U+K is the sum of the potential and kinetic 
energies of the system; p, q are the momenta and co- 
ordinates of the particles. We expand the free energy 
density T/V up to a constant, again with respect to 
the actual configuration, in terms of this strain tensor 



where 
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In the athermal limit T — >■ the free energy J- re- 
duces to the potential energy U , then the expansion 
(m reads 
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Taking the limit T — > of free energy derivatives 
is equivalent to taking the same derivatives, but of 
the potential energy, while satisfying the constraints 
/ = —VU = [10] ■ In other words, the meaning of 
constrained derivatives is that variations with e keep 
the net forces at zero, i.e. 
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such that f EE~VU = 
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We denote these constrained derivative as -^L, then 
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where derivatives in the above expression should be 
understood as taken at e = 0. Satisfaction of the 
constrained derivatives is achieved by allowing for an 
additional displacement field X to the imposed one, 
which we refer to as the non-affine displacement field. 
So, upon imposing a deformation via some transfor- 
mation H, the total variation in coordinates can be 
split into a sum of two terms 



(10) 



where the 's are added to the imposed deformation 
to satisfy the constraint ([5]). The physical meaning 
of this procedure is that the derivatives are always 
computed after minimizing the energy. The origi- 
nal deformation brings the system away from a lo- 
cal minimum on the potential energy surface and the 
minimization amounts to an additional non-affine dis- 
placement which must be taken into account. Note 
that the non-affine field should be found explicitly, 
in terms of potential energy derivatives, and this is a 
part of the derivation below. 

The presence of constrained derivatives calls for 
some care in the calculation. For the benefit of a 
novice reader we present in Appendix |A] some intro- 
ductory remarks to facilitate further reading. 



B. Non-afHne Velocities 

The non-affine corrections X of Eq. (fTO|) play the 
role of y in equations (|AHIA4[) . as their evolution is 
dictated by the constraint ([5]) . The role of x in equa- 
tions (|AHIA4p is played by the strain e, such that in 
complete analogy with (jA4l) . the constrained partial 
derivatives with respect to strain can be written as 
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where partial derivatives with respect to strain should 
be understood as taken at constant X, and par- 
tial derivatives with respect to coordinates should 
be understood as taken at constant e. The quan- 

titles ^ can be calculated by demanding that 
the constrained derivative of the forces vanishes, i.e. 



9fL 



0. We use ([TTt and write 



deaf 



dfl 
deaB 



dXl 



de 



a/3 



dfl 
dXl 



(12) 



We identify the first term on the RHS of the above re- 
lation as the (negative of the) mismatch forces S 



KaP 



jru 



jg—r, which can be calculated explicitly 
from the potential energy, see Appendix [O for the 
case of pairwise potentials. Notice that derivatives 
at constant e can be equivalently taken with respect 
to the coordinates x^^ or the non-affine displacements 
X^, following Eq. ([TOl) . The mismatch forces arise 
due to the imposed affine deformation, before the ac- 
tion of non-affine correcting displacements. The sec- 
ond term on the RHS of (fT2)) contains the derivative 

§xi = ~dxBxl = ""^^^^t- ^^^'^^ negative of 

the Hessian. Using these definitions of the mismatch 
forces and the Hessian, Eq. (IT2]) can we rewritten as 
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We now define the non-affine velocities V] 



OX' 



, which can be calculated by inverting Eq. p^ : 
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Notice that the Hessian is not generally invertible due 
to zero-modes that arise from the translational invari- 
ance of the potential energy; to overcome this we ex- 
pand the RHS of (fT4|) in eigenfunctions of the Hessian 



Finally, 
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Notice that the first order elastic constants do not 
contain a relaxation term; Noticing the stress er — Ci, 
this means that the non afRne displacement field does 
not relax the stress. For an explicit expression of the 
potential energy derivatives in the case of pairwise 
potentials, see Appendix [Cl 



Ka^J = - > . ^ 
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D. Second order elastic constants 



where {ipi)l is the v component of the j'th particle 
contribution to the £'th eigenfunction, Xg is the corre- 
sponding eigenvalue, and the sum over £ should not 
include the zero-modes. 

With the definition of the non-affine velocities, the 
constrained partial derivative (1111) reads 



d 



which is the form that will be used from this point 
on. 



The second order athermal elastic constants C. 



are given by 
we write 
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, evaluated at e = 0; first. 
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where the second equality follows from (1^(71) . We use 
the rule (flGl) and calculate 



First order elastic constants 



The first order athermal elastic constants are 
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Equation (|10l) implies that forces can be calculated 
by taking derivatives with respect to either set of co- 
ordinates, i.e. 
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which can be used in calculating the constrained 
derivative 
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The quantities c, ^ a — and 2* = ^via — can be 
directly calculated from the potential, see Appendix [Cl 
for the case of pairwise potentials. Finally 
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The second term on the RHS of the above equation 
vanishes due to the constraint / = as appears in 
and we are left with 
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E. Third order elastic constants 



The third order elastic constants C. 
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are given 
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We notice now that according to (|16p 
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Also 



— n — = a""'^ = a""" with thesc relations, (l24l) becomes 
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The above expression requires knowledge of the quantity 



for which an explicit expression is unavailable, 



but can be obtained by demanding — ^ 
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Now, 
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Be 



VI 



BE 



dXl 



and — - — - 



<^KX 



5^ , 

Be,^ + ^""^^ BXl ■ 



BE 



de, 



+ vi 



BV: 



rjaP 



With the identity 
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, we invert (|28| to get an expression for 
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Inserting this relation back in ((26|). together with ([141) '^^ arrive at the expression for the athermal third order 
elastic constants 
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III. THERMAL ELASTIC CONSTANTS found in the literature, see for example jlOf for the 



The derivation of the thermal elastic constants, 
in contrast with their athermal counterpart, can be 
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first and second order objects. We provide below also given by derivatives of the free energy with respect 
the third order constants. to strain, see Eq. [5] We first expand the total energy 

The elastic constants at finite temperatures are up to third order in the strain 
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Explicit expressions for kinetic energy derivatives can be found in Appendix [B] and potential energy derivatives 
for the case of pairwise potentials are available in Appendix [CJ We expand the Boltzniann factor 
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The partition function can now be written as 
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With these definitions, we extract the following relations 
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where triangular brackets denote equilibrium averaging. 



A. First order thermal elastic constants 

The first order free energy derivative with respect 
to strain is 
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With the expansion and relations this 

derivative at e = is 
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The first order elastic constants are thus given by 

Cf^^iEf). (39) 



With the expansion ([34]) and relations (p6l) . this 
derivative at e = is 
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where AA = A — (A) for any quantity A. The second 
order elastic constants are thus given by 
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B. Second order thermal elastic constants 



The second order free energy derivative with re- 
spect to strain is 
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C. Third order thermal elastic constants 



The third order free energy derivative with respect 
to strain is 
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With the expansion ([M)) and relations this derivative at e = is 
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The third order elastic constants are thus given by 
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IV. T ^ LIMIT OF THERMAL ELASTIC 
COEFFICIENTS 

We first work out a general expression for the low 
temperature expansion of the mean of any quantity 
A{x) which depends only on coordinates; we start 
with the definition 
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where Zc = J e~T dq. With the notation 5x\, = xl — 
x]^, we expand A around some local minimum x, up 
to fourth order in coordinates 
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Wc now expand the potential energy up to fifth or- 
der in coordinates around the local minimum x, with 
7", Al, Q denoting the third, fourth and fifth order 
derivatives of the potential energy with respect to co- 
ordinates, respectively, evaluated at x: 
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for every x\,. At low temperatures the virial theorem 
insures that the second order term in the potential en- 
ergy expansion is proportional to T, hence the third, 
fourth and fifth order terms are of higher order in tem- 
perature; thus, in the athermal limit, we can expand 
the corresponding parts of the Boltzmann factor: 
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We have omitted the other sixth order term since it does not eventually contribute to the temperature expansion 
of (A) . Inserting the above expansion and (|47)) into (l46l) , we carry out the Gaussian integrals to obtain the low 
temperature approximation of the equilibrium mean of A up to terms of 0{T^): 
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where Ax, Axx etc. denote the tensors of first, second etc. derivatives of A with respect to x. Using this 
expression, we derive relations for fluctuations. Assume that also B and C are functions of coordinates x; with 
the definition A A = A — (A) , we have 
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and 
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In this section we will make use of the following definitions, (see also ([32]) ): 
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A. First order elastic constants — athermal limit The second term in (1571) is 



The first order free energy derivative with respect 
to strain is given by (see Eqs. ([32 |) . (|38l) ) 
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In the athermal limit the kinetic term (jB8|) vanishes 
and following relation ([50]) we are left with 
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From here the athermal limit of the first order elastic 
constants is 
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in agreement with (PT|) . 



B. Second order elastic constants — athermal 
limit 



The second order free energy derivative with re- 
spect to strain is given by (see Eqs. (|32t . (|4T|) ) 
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the athermal limit the kinetic term (jB9|) vanishes, and 
following ([50]) we are left with 
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The second term of the RHS of the above relation is 
proportional to T^, so we are left with 
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C. Third order elastic constants — athermal 
limit 



The third order free energy derivative with respect 
to strain is given by (see Eqs. ([5^ . pil) ") 



T T T 

The first term on the RHS of dH]) is {E^^"'^'^^) = {U^'^"''''^) + {K^"'^''^)] in the athermal hmit the kinetic 
term (|B10I) vanishes, and fohowing (15(11) we are left with 
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since similarly to the second order athermal limit case, it is easy to verify that all mixed terms envolving products 
of kinetic energy and potential energy derivative cancel. The triple product {AK'^^ AK^"^ AK^^) ^ T^, so we 
are left with 
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Finally, we use relation (1511) for the remaining terms: 
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Combining results (|62l) .(j64 l) . (|65p . we arrive at the final result 
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in agreement with (|5ni) . 



V. EXAMPLES OF APPLICATIONS of important issues regarding elasto-plastic behavior 



To justify and motivate the calculation of the non- 
linear elastic constants we present now two examples 
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FIG. 1: Color online: Stress-strain curves. Left panel: 
starting the experiment from a freshly prepared sample 
results in a symmetric trajectory for 7 —7. Right panel: 
starting the experiment from the zero-stress state with 
7 = 70 results in an asymmetric trajectory, see text for 
details. Data was averaged over 500 independent stress- 
strain curves at T = 0.01 where temperature is measured 
in units of e/fcs, see ^ for details 
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0.1 0.2 0.3 0.4 0.5 

7 



FIG. 2: Color Online: Upper panel: Trajectories of stress 
vs. strain for four different temperature at the same strain 
rate 7 — 10"'*. Lower panel: the corresponding values of 
B2 as a function of strain. Data was averaged over 1000 in- 
dependent stress-strain curves at each temperature. Note 
that B2 is negative even when the averaged stress-strain 
curve has a positive curvature, see text for discussion. 



in amorphous solids that cannot be discussed with- 
out invoking these nonlinear constants. The first is 
plasticity-induced anisotropy and the second is the 
predictions of plastic failure. 



A. Plasticity-induced anisotropy 

A freshly produced amorphous solid is isotropic, 
and as such presents a symmetric stress vs. strain 
curve for positive or negative strain. This is not the 
case for the same amorphous solid after it had been 
already strained such that its stress exceeded its yield- 
stress where plastic deformations become numerous. 
This is demonstrated in Fig. [T] A typical averaged 
stress-strain curve for a model amorphous solid (see 
Ref. [4*1 for numerical details) starting from an ensem- 
ble of freshly prepared homogenous states is shown in 
the left panel, with a symmetric trajectory for posi- 
tive or negative shear strain. Once in the steady flow 
state, each system in the ensemble is brought back to 
a zero-stress state, which serves as the starting point 
for a second experiment in which a positive and nega- 
tive strain is put on the system as shown in the right 
panel of Fig. [1] Even though the initial ensemble is 
prepared to have zero mean stress, the average trajec- 
tory is now asymmetric. This phenomenon, sometime 
referred to as the Bauschinger effect |Tlj, shows that 
the starting point 70 for the second experiment (re- 
ferred below as the Bauschinger point) retains a mem- 
ory of the loading history in some form of anisotropy. 

To shed light on the anisotropy of the Bauschinger 



point we choose to measure the sum 



^2(7*) = lim 



(j^VVV 



^xyxyxy 



= lim 



(67) 

which can be determined using the results Eqs. ([^ 
and p5l) . In particular we note that B2 is identically 
zero in an isotropic ensemble. It gains a nonzero value 
when plastic events take place and begin to build 
anisotropy. In Fig. [5] we present results of numeri- 
cal simulations in a typical model of an amorphous 
solid (for details cf. [J|), and present the measured 
value of B2 along the trajectory shown in Fig. [T] In 
addition to the very low temperature trajectory in 
Fig. [1] we show also measurements of B2 for simu- 
lations performed at other three different tempera- 
tures. For all these trajectories B2 was measured 
by quenching to zero temperature configurations sam- 
pled along the thermal trajectory. The magnitude of 
the Bauschinger effect goes down when temperature 
goes up, and this is in good agreement with the value 
of B2 at the Bauschinger point which also decreases 
when temperature increases. We thus see that the 
values of the nonlinear elastic constants can serve as 
natural measures for the degree of anisotropy that is 
built up in an amorphous solid due to plastic deforma- 
tions. Further discussion of this measure can be found 
in Q . Note that in [4j] the value of B2 was obtained 
directly from stress vs. strain curves, computing the 
derivatives numerically. With the results obtained in 
this paper we can compute B2 or any other elastic 
constant directly from the particle interactions. This 
will open up in the future a possibility to define local 
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values of the elastic constants, providing maps of B2 
or other nonlinear elastic constants. 



B. Predicting plasticity 

Imagine an amorphous solid under a given state of 
strain. Can one predict how much additional strain 
is needed to reach plastic failure? Recently we have 
shown, cf. [H, that an accurate predictor of plastic 
failure in an amorphous solid can be constructed with 
the help of the higher order derivatives of the po- 
tential function. Our findings not only offer a pre- 
dictive tool for the onset of failure, but also point 
out the importance of nonlinearities, and in particu- 
lar those that couple nonlinearly "softening" regions 
with strain at larger scales. It is indeed this non-linear 
interaction which produces a reinforcing mechanism 
leading ultimately to a catastrophic event in the form 
of plastic yielding. 

To fix ideas, imagine a simple shear deformation 
applied to a given piece of amorphous solid (for sim- 
plicity in 2D, with immediate extensions to 3D). A 
small strain increment corresponds to a change of 
the i'ih particle positions 



Xi H 
Vi , 



In athermal quasi-static conditions {T 0, 7 — >■ 
0), the system lives in local minima, and follows 
strain-induced changes of the potential energy sur- 
face. Therefore, the particles do not follow homoge- 
neously the macroscopic strain, and their positions 
change cls X-I — y X ' -|- Xi, where Xi denotes non-affine 
displacements. Around some stable reference state at 
7 = 7*, the field Xi, the system energy, and internal 
stress axy are smooth functions of 7. We can thus 
write: 



.(7)^ 



00 



7 , i>„ = Imi — z — - 

' ' ' T->0 d7" 



As the strain increases, the system must eventually 
lose mechanical stability; the "elastic branch" on the 
stress curve ends in a discontinuity as the system fails 
via a first subsequent "plastic event". It is precisely 
at this instability, say at 7 = 7p, that the function 
a^y (7) loses its analyticity. Accordingly we recognize 
that the radius of convergence of the series ((68|) is 
precisely |7p — 7o|, where 7p can be larger or smaller 
than 7*. 

Mechanical instabilities are associated with the 
vanishing of an eigenvalue of the Hessian, which we 
will denote as Ap. Here, as opposed to the analysis 
of 0, Q , we do not make any assumptions about the 
functional form of Ap(7). The vanishing of Ap implies 
that any terms of the elastic constants that are com- 
posed of the inverse of the Hessian Ti^^ will diverge 



at Ap; hence, the strongest diverging terms are those 
which are composed of the largest number of 'H~^'s. 
Accordingly, higher order elastic constants will con- 
sist of stronger divergences, as can be understood, for 
example, by comparing Eq. ([23]) with Eq. (pOl) . Near 
7p the most diverging terms will dominate over the 
rest, so we will consider in the following analysis the 
most diverging terms of each order of the athermal 
elastic constants. 

We begin with Bi = C2^^^, which is given by 
Eq. (1231); following the discussion above, the most 
(and only, for C2) diverging term is Vg^yEl^y/V (re- 
call that V consists of see Eqs. ([HI) , ([15]) ) . Close 
to 7p, the diverging term will dominate, so we can 
write 



Bi 



1 

a; 



(69) 



We continue with B2; the most diverging term in B2 

^3 

this is 



is the most diverging term in C^^^y^y-^ from Eq. pOI) . 



B2 ^ '^l'r,yi'r,yn^i,Tft} 



1 



Bxy'' (^xy'^ pxy I eCp \3 
Ap 



(70) 



Since B2 = we obtain a differential equation for 



Ap(7): 



1 d I \ 

X% c?7 \Ap 



1 dXp 



We re-write this as 



dXp 
c?7 



1 

A^ ' 



(71) 



for which, together with the boundary condition 
Ap|-yj5 = 0, the solution is 



Ap %/7p ~ 7 : 



(72) 



in agreement with 0, [H| • 

With this result, we are able to derive expressions 
for the diverging terms of B3 and i?4, as a function of 
7p; starting from Bi = we obtain the relations 

3a 



B. 



— and B4^ 



15a 



4(7p-7)^ 8(7p-7) 
Solving for 7p, we obtain the prediction 

5B3 



-• (73) 



7p = 7 + 



2B4 



(74) 



where the most diverging terms of B^ and B4 should 
be considered. Notice that one could, in principle 
derive expressions for 7p involving lower order elastic 
constants; see [1] for discussion. 
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We finally derive expressions for the most diverging 
terms of and B4 ; starting from the most diverging 
term in C3, we take another constrained derivative: 



From here 



|-([vvv]-r) 

oe 



= - [T-vv]-n-^ -[T-VV] 

+2 symmetric terms , (75) 



where here and in the following, the contractions are 
only over indices and components of particle coordi- 
nates, and we only consider here the most diverging 

term of ^ . see Eg. ([29]) . 
/ 

Taking another constrained derivative of (|75t re- 



quires an expression for 



which can be 



obtained by applying the rule on 71 ^ 



-in 



a0 



de 



^TaP ' TKX 



-l\mj 
'XV 



(76) 

With this relation and (|29l) , we carry out another con- 
strained derivative of (|75p . keeping only the most di- 
verging terms 



^{[T-vv]-n-' -[T-VV]) 



= r ■ [{n^^ ■ [T ■ VV]) V [n^^ ■ [T ■ VV])] + 14 symmetric terms 



This leads us to the final expression for the instability strain 7p 

,-1 



1 

7P = 70 - 7: 



2 r • [{n-' ■ [T ■ VxyVxy]) Vxy ■ [T ■ VxyVxy])] ' 
I 



(77) 



Notice that the above expression for 7p can be cal- 
culated numerically by solving two linear equations; 
first, for V using With the solution for 

hand, one can then solve equation (|28l) for 

'H~^ ■ [T • VV). How these predictions work in prac- 
tice can be read in Q. 



sion of the evaluation of these objects to local coarse- 
grained fields. Such an extension will be presented 
in a forthcoming article. We demonstrated the use 
of these nonlinear objects in the context of under- 
standing the plasticity-induced anisotropy that arises 
in amorphous solids after a straining trajectory and 
in predicting the plastic failure of amorphous solids 
to increasing strain. 



VI. SUMMARY AND CONCLUSIONS 



In this paper we derived closed-form expressions for 
the nonlinear elastic constants of amorphous solids up 
to third order. We presented both the thermal and 
the athermal theory, and demonstrated that the lat- 
ter is obtained as a limit of the former when T — >■ 0. 
The expressions derived above should be useful in nu- 
merical simulations where knowledge of these nonlin- 
ear constants is indispensable due to the high val- 
ues of the stresses obtained near mechanical instabil- 
ities. In particular these expressions will allow exten- 
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Appendix A: Constrained Derivatives 



Appendix B: Kinetic energy derivatives 



To provide an intuitive demonstration of the mean- 
ing of constrained derivatives, consider two functions 
of two independent variables h{x,y) and g{x,y). We 
want to understand, for instance, how one should take 
the partial derivative of h with respect to x, keeping g 
constant. The constancy required from g(x, y) means 
that y now depends on x, so that variations of x must 
result in variations in y, namely 



, dy 

dy^^ 
ox 



dx 



We write the total variation of h as 



dh = 



dh 

dx 



dx 



dh 
dy 



dy . 



(Al) 



(A2) 



We now impose the dependence of the variation in y 
on the variation in x: 



ox 



dx + — 



dh 
dy 



dy_ 

dx 



dx 



From here 



dh 




dh 




dh 


dy_ 


dx 


g 


dx 


H 

V 


dy 


X dx 



(A3) 



(A4) 



It is important to notice that in this framework, con- 
strained partial derivatives do not necessarily com- 
mute; to exemplify this issue, consider the constraint 
g{x, y) = x"^ -\- y^ = 1 m. the range of positive x and 
y. The partial derivative of y with respect to x at 
constant g is 



dy_ 

dx 



We first consider 



dx"^ 



d_ 

dx 



Next, it is immediate that 
dx"^ 



d 


dy_ 


v,g ~ 


g 



d_ 

dx 



= 



(A5) 



1 

V ' 
(A6) 



(A7) 



since obviously any variation of y keeping y constant 
is zero. So, generally 



dx'^ 



v,g 



^ dx^ 



We finally derive expressions for kinetic energy 
derivatives that appear in ([32|). Given a transfor- 
mation of coordinates H , the transformation of mo- 
menta is dictated by requiring that the transforma- 
tion is canonical, namely that 



dx'^ dp'p dx'^ dp'p 
dxv dpu dpi, dxi. 



■JalS 



(Bl) 



where prime denotes transformed coordinates and 
momenta, and in this chapter repeated indices are 
summed over. Assume that = Ap^Pn, such that 

1^ — Ap^S^i^. Inserting this in equation (jBl|) : 

Ha^S^^Ap^S,^^ = HavAfj^ — HauA^p = 6ap ■ (B2) 

This means that HA^ = I, or A — (iJ^^)^, such 
that the transformation of the momenta is 

Pa - {H'^ZpPP ■ (B3) 
The kinetic energy after imposing a small strain is 

K = \pipl^\{H-')l,pl{H-T^pfp 



hH^IiH-')lpPlp}3 



(B4) 



From ([U we have H^H = 2e + I, hence 

H-\H-^f = {H^H)-^ = {2€ + I)-^ . (B5) 

Expanding this up to e'^: 

(2e + I)-^ ~ / - 2e + - 8e^ + 0{e^) . (B6) 

Inserting this back to ([B4|: 

K ~ ^SaflPaPp - eaPPaPp + 2eau£^fiPaPp 

(B7j 

- 4:eai,ei,rjerjPPaPp + 0{e ) . 

We can now calculate 
dK 



dea 



'iep^e^rjPliPa - '^^au^PvPlp''^ - '^^uri^qaPlp'p 



(B8) 



(A8) The second (symmetrized) derivative is 



g,v 
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- 2Sf3„e.n0plpl - 2e/3;.pJ,K - 2Si3r,(-uep},pl - "^^PvPlPa (B9) 

- '^Sav'^nep'^pl - 26p^ea0plplj - 2(5„,,e^0pJ,pg - 25p^eaePePl 

- 2ejjaPlP0 - "^^aveeuPeP'is ~ '^<^uaP\,p'i3 - '^Sa,,eefjPeP'i3 , 



and the third (symmetrized) derivative is 



- Sl3„SKr,p]^Pa 


- 6^05xvp\^p'„ 


- 5pr,5KvP\Pa 


- 6^0dxr,plPa 


- S/sJxvPlPa 


- Sx^S^^pIipI, 


- SpnSxi^PlPa 


- S^pS^nPlPa 


- SauSKl3Pl,Px 


- SlSuSnaPxPl^ 


- 6a-q5f,l3Plp\ 


- Sl3r,SKaPxPl 


- 5a^5^i3pl,pi 






- 60r,SxaPipi 


- S^nSxaPlp} 


- 5a'n6xuP\p''f, 


- 6f,y5xaPl,p''0 


- Sai^Sxr^p'^p} 


— S-^n^naPvP'/B 


~ ^ari^KvP'xP'^ 


— ^xv^KaP\p^p 


~ ^av^KT^P^xP^P 



(BIO) 



Appendix C: potential energy derivatives for 
pairwise potentials 

In this appendix we calculate all partial derivatives 
of the potential energy with respect to strain and par- 
ticle positions, of all required orders. In the following 
section, partial derivatives with respect to coordinates 
X should be understood as taken at constant e, and 
partial derivatives with respect to strain e should be 
understood as taken at constant non-affine fields X 
(no relaxation allowed for). 

We first carry out partial derivatives of the poten- 
tial energy with respect to particle coordinates. The 
potential energy is given by 



(CI) 



where (j)^^ is the pairwise interaction potential be- 
tween the i'th and j'th particles. The negative of 



the forces are given by 



J a 



dxi 



dr^^ dxi, 



(C2) 



9x1 



Since r*-' = ^ (a;^ — a;^)(a;^ — a;^) then 

jTj-{5^^ — (5*^). With the notations (pr,4>rr and (j)rrr 
for the first, second and third derivatives of 4>{r) with 
respect to r, respectively, we now have 



£f = ^ E ^(^^' - ^'') = 5 E ^^^^ - 1 E '^^^■^ = ^ E ^^^^ + 5 E ^^/^ = E (C3) 

ijte e^j i^t m i^e 

The Hessian is defined as 

'^^ dx^^dxl, 2 dx^^ dr^^ dxl^ 2 dx^^ i^, r"^ " ^ ' dx^^ ^ r« « 

= E (t^ - rfrfi^'^ - 5'^) + S., E ^{S- - S'^n ■ (C4) 

The off diagonal term i ^ j of the Hessian is given by 

njij _ { Vrr Vr \ tj i] X /prN 



15 



and the is given by 



Q/3 



(C6) 



We now calculate 



assuming first i ^ we take a derivative of Eq. (jC5 



0rrr ^^rr , \ ij ij ij . ( 0rr \( s: i A _l A 

(-^»i)3 " (-^»i)4 + (-^»j)5 y ''/3 ''X (-^»i)2 " (-^»i)3 ^ +f'/3x?'a +0a/3r'x 



(C7) 



If i ^ j =^ £ then 



^ ^ ^ . = 0, so in the above expression £ must be equal to either i or (but not to the 



d-'U 



both). Due to the symmetry 
first, if two of the particle indices are equal, and different from the third, then 



dx^dx-j^dx^^ 



dxl^dx^dx-'^ 



dx-^^dxl^dxl 



-, we can limit the discussion to two cases; 



(j-iiY (r'J)'* (j-iiY 
where we assumed i ^ j- The diagonal term is then 



dxi^dx'pdx). 



(C8) 



(C9) 



We now turn to derivatives with respect to strain e; consider first the change in distance 6r'^^ between 
particles i and j, before and after imposing a deformation. For a given imposed deformation represented by 
i?, each component of the pairwise distance transforms via — >■ Hapr^^, so the distance between the pair 

(«, j) after the deformation is H^^^Hypri r'^ . Since by definition H^H = 2e + 1, the change in the pairwise 
distance is 



~ r 'J ( 1 + 



Tp r 



(CIO) 



r^'^3 ip^^ 



+ 2- 



We write the potential energy as a sum of pairwise contributions U = X]i< j j ^'^'^ expand the pairwise 
potential in terms of Sr"^^ : 

= + d^^Sr'^ + \(lfl5r'3 5r'^ + \(lfl^5r'^ Sr'^ 6r'^ + . . . , (Cll) 

Plugging (|C10p in the above expansion gives (j)^^ in terms of e, keeping terms up to e'^: 



I ±^ c fr'^3 rp'^3 ,Y''^3 ^"^3 <Y''^3 

-r Q*=^ctf3^iy7]^Kx' a ' j3 ' ' r] ' k ' x 



(j.ij^5 (r*J)4 [r^jy 



(C12) 



From here the derivatives of (j) with respect to strain can be calculated as 
de^p~ r-i ^Wa'-^^.^ (^^^,,^2 



If f r^j r^l r^j r^i r"^] r"^] 
2 '^vri'^Kx' a ' p ' V ' r) ' K ' X 



(r«-J )3 J ' 

(C13) 
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dey^deap 



and 



f: j-ij r''] r^j r^J 



(jij^ (7'ii)4 [rijy^ 



— r'^j r'^i 



f 3<j)l? Z<j)% <j)%^ 



The derivatives with respect to strain of the total potential energy at e = are given by: 



dU 



de. 



and 



E 



52 iJ 



i<3 



13 r^3 „t3 r^3 



I r'rr 



denxd^v-qdecxp 



P 7] K X 



i<j 



mi 



Turning now to the mixed derivatives, we first consider g^kg^ — ! following (jC13|) . we obtain 



92 



Next we calculate the mixed derivative 



[ ^rr _ Vr \ 13 i] i3 j.ZJ_(S. ^-'J 1 A ar^l \ 



following (jC14p . we obtian 



( 



mr 



j,^3 ^tj „U 

X P 1 



(§3k _ 

\ (r*-')3 (r*^')'* (ry)'5 
/ ■ ■ ■■\ / 6*^ 6""^ 

X r^3 i3ri3 iX., r^3 r^3 r'^3 _i_ X r^3 ^^3 13 _i_ X ^^3 ^^3 13 ] ^rr Yr 

yax'p'iy'ri ^"Px'a'u'ri ^ "vx' a ' p ' 1] ^ "vx' a ' p ' u J \^ij^2 (^j.ijy 



The last mixed derivative is 



following Eq. (IC18I) . we obtain 



dxi^dx^deafz 



r^3 r^3 „i3 „t3 
rp ^x 



+ 



H JT {SuaSpx + SupSax) 

T 



For any quantity A^^ such that A*-^ — — ^-^^ we have 



ik 



(C14) 



(C15) 



(C16) 



(C17) 



(C18) 



(C19) 



Y^{^53k _ pk^Ji^J ^ 1 ^ ^(JJ/^ _ 5''^)A'^ = 5 E " ^ E ^''^^ E ^ 

Since in equations (jClSp . (|C19P and (|C20p the terms that correspond to A^^ in the above equation are anti- 
symmetric in so we can directly write the results for the mixed energy derivatives: 



dx^Be. 



aP 



-E 

e=0 j^fc 



J,ik lik \ lik 

r'rr T^r \ ik ik ik , r^r ( ^ ik , e ik\ 

pfejl ~ ■pfcjl ) rp + —k['^'^arp + Oi^pr^ ) 



(C21) 



dx^de^ridcap 



E 

i=ik 



30: 



X. r^3 „ij „ij 
^ax' p ' V ' 7] 



ik 



^j,ik'^3 ^j,ik'^4 ^j,ik'j5 



^tk fy,ik ^ik ^ik ^ik 
^X a ^p ^ri 



X,. r^3 r^3 r^3 \ X r^3 ^''^ ^^3 \ X r^3 r^3 „ij 



l^ik 



lik 



^j,ik'^2 ^j.ik'^3 



(C22) 
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and 



dx^dxldeap 



4,ik jAk 



<i3f„ 30f^ 30f 



„ik ih ih ik 



In the case for which t ^ k, this is 
93 f/ 



dx'idx^deap 



rrr 



30rr , 30^*^ \ fk„ek„ek„ek 



(^j,lkY (^j-tk^ (j-lky 

ilk j,£k 



/fc/fc4.A /fer-^fe4.A /'^/fc+A 0/'=/'=^ fC24) 



and the diagonal part is just 



dx^dxldeaf: 



e=o k dx^x^x^de^fi 



(C25) 
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